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Abstract 



o 

• Multi-frequency interferometry (MFI) is well known as an accurate phase-based measurement scheme. 

■ The paper reveals the close relationship between the unambiguous measurement range (UMR), the outher 
\ and the MSE with the frequency pattern in MFI system. We note that the theoretical rigorous UMR of 

MFI deduced in the literature is usually optimistic for practical application and derive a more practical 
expression for the UMR. It is found that the least-square (LS) estimator of MFI has a distinguished 
^ ■ "double threshold effect". Distinct difference is observed for the MSE in moderate and high SNR region 

(denoted by MMSE and HMSE respectively). Beside the conventional threshold effect caused by outher 
. in low-SNR region, another threshold effect occurs during the rapid transition from MMSE to HMSE 

■ with increasing SNR. The closed-form expressions for the MMSE, HMSE and CRB are further derived, 

(N 

■"sj" ■ with HMSE coinciding with CRB. Since the HMSE (or CRB) is insensitive to frequency pattern, we 

_ focused on MMSE minimization by proper frequency optimization. We show that a prime-based frequency 

. interval can be exploited for the purpose of both outlier suppression and UMR extension and design a 

(N 

special optimal rearrangement for any set of frequency interval, in the sense of MMSE minimization. 
An extremely simple frequency design method is finally developed. Simulation and field experiment 

■ verified that the proposed scheme considerably outperforms the existing method in UMR as well as 
, MSE performance, especially in the transition from MMSE to HMSE, for Gaussian and non-Gaussian 

channel. 
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I. Introduction 

Multi-frequency interferometry(MFI) is an accurate phase-based ranging method, widely used in high- 
accuracy ranging m |[23l . localization IJl or deformation and ground displacement detection ll33l . which 
estimates the range using the phase measurements recorded at multiple frequencies. 

How to select the measurement frequency in MFI system? There mainly exists three kinds of criterion. 

• Maximizing the unambiguous measurement range (UMR). 

• Minimizing the mean squared error (MSB) of range estimation. 

• Decreasing the outlier(the wrong estimate lies outside the main lobe of the cost function [17] ) 
probability and enhancing the robustness to noise. 

The paper will reveal the close relationship between these three criterions with the frequency pattern 
(refer to the frequency spacing of adjacent frequency in the paper) and pay special attention to the 
optimization of frequency pattern. 

A. Related work 

Equally spaced frequency is employed in RIPS |12|1 and the UMR is the synthetic wavelength of adjacent 
frequency, i.e. c/Af. In the geometric series of wavelengths is used and the UMR is extended to 
the synthetic wavelength of the closest two frequencies or the largest synthetic wavelength. |[23]| adopts 
the beat wavelength coincidence method and the UMR is further enlarged by a factor of positive integer 
relative to that of H. 

Ranging accuracy of multi-frequency interferometry system is another important target to be optimized. 
However, there exists limited works concerned with the accuracy improvement by frequency optimization 
in the literature. 

Since the estimation accuracy of RIPS method is in proportion to frequency separation, while the 
UMR is in inverse proportion to frequency separation HI m. Then, there exists a compromise between 
the accuracy and UMR. The multistage beat wavelength method is exploited for range estimation and 
the measurement phase noise will be amplified by Aoi/Ao,j+i when two neighbor synthetic wavelengths 
Aoi, Ao.i+i are used for phase unwrapping(determine the integer number of wavelength) in each stage 
in 1231 . In order to maximize noise immunity, the wavelength ratio should be equal and the phase noise 
will be uniformly amplified in each stage. With this in mind, a geometric series of wavelengths is formed 
m. However, the optimization criterion of |4| is built on the beat wavelength method. Although it is a 
fast estimation method, it is far from the optimal one in term of estimation accuracy since only partial 
information is used in each stage. 
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Note that similar problem arises in array design, that is, how to improve DOA estimation performance 
by optimizing the antenna positions of the n-element linear antenna array. The linear array geome- 
try(antenna separation) design is most similar to our frequency pattern(frequency separation) design. The 
antenna positions are optimized to minimize the CRB of DOA estimation by numerical optimization in 
|[3l . Based on the CRB principle, ifTTI also considers the linear array optimization problem for joint range 
and DOA estimation. 

However, it seems not to be reasonable to optimize array geometry or frequency pattern only relying 
on CRB since it is a local measure of estimation accuracy. The ambiguity function is introduced when 
the global accuracy performance and the outlier are evaluated IfTTI ifTll |fT6l ~ |fT9l . 

Nonuniform linear array optimization is discussed in |[T2l with emphasis on the comparison between 
uniform array, nonredundant array and minimum redundant array. Array geometry is optimized in ifTSl 
to minimize the outlier probability. For each given geometry, the ambiguity function reaches its minimal 
at a certain direction and the optimal array is the one with the minimal value maximized. Both CRB and 
the ambiguity function are used to optimize the two-dimension array in |[T9l . and the genetic algorithm is 
adopted to search for the best array with the minimal CRB, under the constraint that the outlier probability 
is below a certain threshold. 

Outlier or sidelobe suppression is also an important topic in step frequency radar. Due to extremely large 
bandwidth, the estimation accuracy is usually satisfied and more attention is paid to outlier suppression. 
To this end, random frequency step or different number of the frequency repetitions are proposed |[30ll 
|[3T1 |[32l . This is quite different from narrow band ranging. It is interesting to find in later section that 
the frequency pattern of our method is quite different from that in step frequency radar, where more 
weight is given to the frequencies distributed in the center of band. 

B. Our work 

It is known that the theoretical rigorous UMR of MFI is the least common multiple (LCM) of all 
wavelengths when expressed as integers by quantization Q lfT4ll llTl . We prove that the UMR is the 
LCM of all synthetic wavelength of adjacent frequency or in inverse proportion to the greatest common 
divisor (GCD) of all frequency spacing when the initial frequency is properly chosen. Thus reducing 
the GCD will readily extend the measurement range. Compared with the wavelength-based method, the 
frequency spacing-based method can make full use of the band in the sense that there have far more 
frequencies to be picked up for UMR extension for a given bandwidth |29|. What is more, with the 
initial frequency fixed, it is possible to adjust the order of the frequency spacing (known as permutation) 
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to generate various set of measurement frequencies maintaining the same UMR since the GCD of all 
frequency spacing remains unchanged. This forms the key foundation of later MSE optimization via the 
permutation of frequency spacing. It is also pointed out that the rigorous UMR is overly optimistic in 
case of random initial frequency and a more practical expression for the UMR is derived. 

Inspired by the method in array design, we also resort to the concept of ambiguity function for analysis 
and alleviation of the outlier probability. Note that obvious differences still lie between the two systems. 
The parameter space of ambiguity function is "bounded" to the direction of —180° 180° for antenna 
array. It is not the case for MFI since the UMR will increase monotonously as the GCD of frequency 
spacing decreases. 

Constrained by the size and cost, the number of antenna is usually limited for array design while a 
large amount of frequencies may be required for accurate ranging in MFI. Therefore, the search based 
method used to cope with the outlier in antenna array could not be directly applied to the MFI system. 
Based on the ambiguity function, we reveal the relation between the frequency pattern and the outlier 
probability and a design criterion for outlier suppression is then provided. 

The relationship between the theoretical attainable ranging accuracy and the frequency pattern of MFI 
is the major issue of the paper. It is well known that the ML estimator exhibits a threshold effect, i.e., 
the MSE increases sharply below a certain signal-to-noise ratio (SNR). Surprising, it is found that the 
LS estimator of MFI system shows a distinguished (Impressive) "double threshold effect", especially for 
fi » B, which is independent of frequency pattern. Distinctive (Obvious) difference is observed for 
the MSE in moderate and high SNR region (denoted by MMSE and HMSE respectively). Beside the 
conventional threshold effect caused by outlier in low-SNR region, another threshold effect occurs during 
the rapid transition from large MMSE to small HMSE with increasing SNR (see figS]). The closed-form 
expression for the MMSE, HMSE and CRB are further obtained with HMSE coinciding with CRB. Under 
the assumption that fi ^> B, which is true in most cases, different frequency design methods have almost 
the same HMSE and it is not the case for MMSE. Therefore, we focused on MMSE minimization by 
proper frequency optimization. 

A quite simple prime-based frequency design method is finally developed with no need of searching. 
That is, if priori knowledge is provided and the range uncertainty can be restricted to the mainlobe region, 
i.e., the outlier is excluded, the optimal measurement frequency is proved to be densely distributed on both 
ends of the band. Otherwise, a prime-based frequency interval is firstly constructed for the purpose of both 
sidelobe suppression and UMR extension, followed by a special optimal rearrangement of the frequency 
interval. The proof of its optimality is given for the first time, in the sense of MMSE minimization. 
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Different from the exhaust search or numerical optimization method, which have prohibitive com- 
putational complexity with increasing antennas or measurement frequencies (Suppose N measurement 
frequencies should be selected from M possible frequencies, there exists combinations of measure- 
ment frequencies, so the antenna number = 4 is used in j TSl ). the proposed method may be free 
of the frequency number constraint. Moreover, for given bandwidth and frequency number, it is usually 
thought that the UMR and the ranging accuracy can not be improved simultaneously 1281 . However, the 
proposed method and subsequent experiments have changed this opinion . 

Apart from the simulation analysis, field experiment in outdoor non-Gaussian channel (multipath error is 
incorporated inevitably) is also performed. The experiment results demonstrate that the proposed scheme 
considerably outperforms the existing method in UMR as well as MSE performance. Meanwhile, the 
simplicity and robustness even to non-Gaussian error are also attractive for its practical use. 

Notation: Upper (lower) bold face letters are used for matrices (column vectors). E[x], l[x], 
Re{x} and Im{x} denote the operation of taking expectation, absolute value, phase angle, real part and 
imaginary part of x, respectively. [x]2tt denotes modulo-27r operation, which reduces x to the interval 
(— vr, vr]. X and x* denotes the estimate and optimal solution of x. Iat denotes the N x N identity matrix; 
In denotes an x 1 all-one column vector; j = ; The mth row and nth column entry of matrix 
A are denoted as A(?ti, n) or [A]^„. The trace of A is given as tr(A) = i')'^' i')^ 

and i')^^ denote matrix transpose, conjugate transpose and inverse operators, respectively. We use script 
letters A to define sets and by |^| its cardinality. 

II. Problem Formulation and Performance Analysis 

A. System model 

In the absence of noise, the measurement phase of multi-frequency interferometry ranging system is 
related to the range by the following equation HI ill HI ll23l 



where Aj = c//j is the wavelength of the carrier frequency /j, and c is the speed of light, go is the true 
range (it is the q-range in radio interferometry and path length difference in optical interferometry)and 
(/7o(i) is the ideal measurement phase wrapped to the principal interval of (— 7r,7r]. The subscript "0" 
denotes the true or ideal value. Equation ([T]) is equivalent to 



(/3o(i) = 27r-- 



(1) 



(2) 
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where is an unknown integer. For single measurement frequency, the range must satisfy l^ol < Ai/2 and 
then qo = ip{i)Xi/2TT with n, = 0. Otherwise, ambiguity appears because rii is undetermined. Therefore, 
multi-frequency ranging is usually needed. 

There are various methods to estimate qo with a set of phase ip{i) given, such as the CRT method [8 | 
|[T3]| lfT4l . the LS-based search method (21, the modulo conversion method |[TOl . the multistage unwrapping 
method llH Q ll2ll (23} . the excess fractions method [24|. Among which, the LS-based search method can 
achieve best estimation performance and will be adopted in the paper (we concentrate on the frequency 
design rather than estimation method), in the form of 

N 

q = argmin S{q), S{q) = J2 iM^) " '^<?(^)]2^)^ 
Q 



i=l 



2tt- 



2n 



27r— -I- I 



(3) 



27r 



where S{q) denotes the cost function and ip{i) is the measurement phase with phase error 6e{i)- 



B. Maximum Unambiguous Range 

Before a formal derivation of the unambiguous range, we begins with an intuitive interpretation of the 
ambiguity problem. 

Fig. [T] shows the relationship between theoretical and measurement phase without noise under different 
range q. The theoretical and measurement phase are represented by red lines and blue stars. Although 
the true range q as well as the theoretical phase are different, the measurement phase are the same with 
A/ = 1.5MHz. Note that the measurement phase is all the information we have in estimating q. In 
other words, the range g = 50 and q = 250 are indistinguishable for A/ = 1.5MHz and the ambiguity 
happens. If we measure the phase with a frequency spacing of A/ = 3MHz instead, the difference will 
immediately appear. This suggests that the ambiguity is closely related to the frequency spacing. 

For multi-frequency interferometry ranging, the phase under q + AL and q are given by 

A = l,2---N 



c 



27r 

Q + AL 



(4) 



2lT 



Ambiguity occurs if and only if the ranges q and q + AL have the same phase in all measurement 
frequencies /j, i = 1, 2, • • • A^. That is. 



= 



C 



(5) 



2tt 
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Then 



AL = niAi = n2A2 



(6) 



where G Z and A,; is the wavelength of frequency /j. Clearly, there are indefinitely many solutions 
of AL satisfying equation system in (H). Obviously, the minimum among those solutions is defined as 
the unambiguous measurement range (UMR). Equation (|6]l implies that the UMR is the least common 
multiple of all wavelengths lIlTI . 

To get closed-form solutions of AL in terms of frequency spacing instead of wavelength, equation 
system in ^ is transformed into an equivalent form as follows: 

• Condition 1: Equal initial phase at initial frequency ipq^^L{l) = ipq{l). 

• Condition 2: Equal phase increment between adjacent frequencies A(pq{i,i + l) = A(^g+AL(^i i + 
where Aipg{i,j) = [ipq{j) - (Pqii)]^^. 

From Condition 1, we get 

q + AL' 



27r/i 



2t7 

AL 



2^/i- 

koc/ fi 



c 

'■ koXi 



2n 
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From Condition 2, we get 



c 



27r L J 27r 

AL = kic/Afi 

where Afi = /j+i — /j (assume /i < /2 < • • • < In) and k^, ki are positive integer. 

Theorem 1: If A/min is the greatest common divisor (GCD) of all the adjacent frequency spacing and 
/i = /cA/min, then the UMR is AL = c/ Af^^^. 

Proof: Without loss of generality, we assume /j+i > /j and the separation of two adjacent frequencies 
is — fi = Afi = fcjA/min, ki is a positive integer. It is clear that AL = c/A/min = he/ Afi, and 
AL = kc/fi, then AL = c/A/mm is one of the ambiguous range. If it is not the minimum ambiguous 
range (UMR), let AL' = AL/m = c/(mA/min) be the minimum ambiguous range, m is positive integer, 
using condition 2, then 

27r 

kiAfrainAL =Pi2TT 

c 

fi+1 - fi = mpiAfmin (7) 

where pi is arbitrary positive integer. ^ indicates that the greatest common divisor is mA/min instead 
of A/min- This conclusion conflicts with the hypothesis. Therefore, we claim that c/A/min is the UMR. 

■ 

Corollary 1: For equal-spaced measurement frequency, i.e., Afi = A/, and fi = kAf, then the UMR 
is AL = c/Af. 

Proposition 1: Suppose A/min is the GCD of all the adjacent frequency spacing, if the condition 
/i = fcA/inin is not met in Theorem [T] and let fi = {ki +e)A/min, —0.5 < e < 0.5, go is the true range, 
then 

N IN 

lim S{q\q = qo + c/Af^^ - Xr' / J] Xr') = 

Proof: See Appendix lAl ■ 
Remark: This impUes that c/A/min — eZ^i^i j Si=i i^ ^'^'^ °f ambiguous range in the 
limit sense. Similar to the proof of Theorem [T] it is easily to verify that it is the minimal ambiguous 
range(UMR). Otherwise, the UMR in Theorem [T] is not the minimum one when e = 0, contradiction. 
If fxjB » 1 and the noise(even extremely small) is considered from the practical point of view, the 
cost function at go + c/A/min — eYld=\ "^^^ j 5I^i=i -^^^ i^ completely indistinguishable from the one 
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at qo or qq + AL, with AL being the theoretical rigorous UMR deduced in In other words, the 
difference between the cost function at qo + c/A/min — ^Y^iLi j Yl!i=\ Qo is only theoretical 

distinguishable. We then define AL = c/A/min — ^J^iLi j ^i=i ^T"^ practical UMR, denoted 

as P-UMR. The UMR in Theorem [T] is a special case of Proposition [T] with e = 0. 

Note that the UMR of CRT method is the product of all wavelengths according to Q when expressed 
as prime number by quantization ||5l S and is much larger than the P-UMR in Proposition [T] it becomes 
an overly optimistic estimation since it is not attainable in practice. 

Since eX^i^i j 5I^i=i ^^'^ < Ai = c//i is negligible compared to cj A^ram under the assumption 
/i » B, the UMR is thought of c/A/min for both /i = kAf^cavL and /i / /cA/min in latter part of the 
paper. 

C. Outlier 

Inspired by the result in |[T8l |[T9l . we introduce the ambiguity function to cope with the outlier(the 
wrong estimate lies outside the main lobe of the cost function) in multi-frequency interferometry ranging. 
It is proved in IITtI |[T9l that the outlier probability is in proportion to the ambiguity function 



where 



C(A<?) 



s(g) 



s((7)^s((7 + Ag) 



iV2 



gj27r/iq/c gj27r/2g/c 



1 

iV2 



N 



(8) 



J2TrfNq/c 



It is obvious that the ambiguity function attains its global maximal 1 when Aq = kAL, where k 
is an integer. This kind of ambiguity(UMR) is inherent and inevitable. The UMR is excluded with the 
constraint |go| < AL. Meanwhile, for < \Aq\ < AL — Bm is the width of mainpeak, there still 
exists the second maximal which may exceed the value achieved at the true go with the help of noise 
and lead to large error, i.e., outlier. The objective is to solve the following optimization problem 



mm < max — 
Af I Aq iV2 



N 



subject to ^<\Aq\<AL-^ 



Since 



1=1 k=i+l 



N-1 N 



(9) 



(10) 
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Then the problem becomes 



I 



N-l N 



mm 

Af I Aq 



maxRe < ^ ^ 



1=1 k=i+l 



Note that we are interested in those Aq satisfying —{fi — fk)A.q ^ 1 , which are prone to resulting 
in outlier. Using a first-order Taylor series expansion of e^, we have 



N-l N 1 N-l N 



i=l k=i+l 

According to Q 

N~l N 

E E 



i=l k=i+l 



— (/,-A)Ag 
c 



2tt 



i=l k=i+l 



2lT 



2tt 



c 



R 



2n 



^A.rAf 

c 



where 



Af = [A/l,A/2,•••A/^_l]^ r 



1 ••• 
1 1 ••• 

: : 1 

1111 



2n 



(11) 



(12) 



R-^ = I 



■Af-l- 



uu 



■, R = I^_i + UU^, U=[l,l,---l]^ = l7V-l 



It becomes the max-mini optimization problem 



— AgAf^r^' 

c 



max < mm 
Af [ Aq 



R 1 



2n 



—AqTAi 

c 



2n 



When is large, it is approximated as 



max < mm 
Af 1^ Aq 



— AqAi^T^' 

c 



2tt 



^AqTAi 

c 



2lT 



(13) 



(14) 



(15) 



We are not intended to find the optimal Af , in term of minimal outlier probability, by exhaust searching 
due to the formidable complexity when the parameter space is large. 
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D. MSE 

Suppose the measurement phase at the i-th frequency is ip{i) = [lirgo/Xi + 0e(*)]27r' estimated 
phase is (p{i) = [27rg/Aj]2^. where qq, q are the true and estimated range, 6e{i) is assumed to be i.i.d 
white Gaussian noise. The LS-based objective function is given by 

N 
i=l 



Since [[x]2n - [y]2n]2n = - yhiT, let Aq = q - qo, then becomes 

n^TTf.Aq 
i=i ^ 



mill > 



27r 



(17) 



c 

We will derive the MSE using the perturbation analysis approach ||26l . The MSE for the moderate and high 
SNR are denoted by MMSE and HMSE in the paper. For SNR high enough, i.e., |Ag| < \m ■ ■ ■ < Xi- 
In this case, ([TT] ) may be simplified as 

Y^IM^.es)] (18) 



mm , , 



The optimal Aq satisfy 



i=l ^ 

2vrE^=i/| 



= TZ^N .0 (19) 



Then the HMSE is readily obtained 

HMSE = E 



For moderate or high SNR, \ Aq\ > \n, but \Aq\ < c/2B, then the error could not be obtained directly. 
Similar to ([TT]) . we may transform the estimator into (since Aq is the one minimizing (ITtI ). then it leads 
to small (27r/jAg/c — 9f,{i)) for moderate or high SNR): 

maxRe j^exp |j - }| (21) 

When fi/B » 1, the problem may be well approximated as (the approximation is reasonable, as is 
confirmed in later simulation, see Fig|2] and Figj4| 

N 



max 



max 
Ag 



2=1 k=i+l 
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85 90 95 100 105 110 115 

range of q (m) 



Fig. 2. Comparison of the actual and approximate cost function 



or equivalently 

N-l N 
y i=lfc=j+l\ L 



(23) 

When \ Aq\ < c/2B, then \Aq\ < \c/2{fi - fk)\. Omitting the influence of 9e{i) - 9e{k) for moderate 
or high SNR, ( [23] ) can be rewritten as 



Define 



c 

*(Ag) = [0A,(2, 1), 0Ag(3, 1), • • • MN, l)f 

@e= [(^e(2) - e,{l)), (0e(3) - 0e(l)), " " " (^e(A^) " 0e(l))]' 



The problem can be expressed as 

N-l N 



min^ (0Aq(i,A:))2=min$^(A(/)R-i*(Ag) (25) 



^ 1=1 k=i+l ^ 



Where R ^, R, u are defined in ([T3] ). The optimal solution obey 

a(*^(Ag)R-i*(Ag) 



5(Ag) ° 
(a [*(Ag)] a(Ag))^ R-i*(Ag) = (26) 
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Since d [0Ag(i, 1)] /d{Aq) = f (/. - /i) = ^ El'Ji ^fk, where A/^ = A+i - f^, and 



27r 

c 

27r 



Af-l 



A/i,A/i + A/2,-- - E 
fe=i 

rAf 



It follows that 



Then 



Since 



APr^ K-'^iAq) = 

Af^r^R-i ^rAf - e. 1 = 



Ag 



E 



Aq' 



E 



c Af^r^R 

2^AFr^R-irAf 



.2 AFr^R-iE 



R-irAf 



47r2 



T 



(Af^r^R-irAf)2 



cj| (Iat-i + uu^) 



a|R 



Therefore, the MMSE for this moderate SNR is given by 

^2 



MMSE = E 



Aq^ 



47r2 Af^r^R-irAf 

^2 



47r2 Af^r'T (I^_i - uu^/N) rAf 



(27) 



(28) 



(29) 



(30) 



It is well known that the phase noise in ([3]) follows the wrapped normal distribution due to modulo 27r 
operation |[27l |[28ll . So the estimation is usually not unbias and the CRB does not exist. But the noise 
can be approximated as normal distribution under the assumption of high SNR. 

Consider the signal y{k) = e-'2'^'?/Afc _|_ ^i^k), where q is the parameter to be estimated and n{k) is the 
complex Gaussian noise with zero-mean and variance E[n2(A;)] = fj2 = cj2 and the signal-to-noise ratio 
(SNR) is defined as SNR = l/cj^. The signal can be expressed as 

y{k) = e^'2-«A'= + n{k) = e^'2-9A.(i + n{k)e-^^'"'/^') (31) 
Let n'{k) = n{k)e^^'^'^'^^^^ = n'^{k) + n'j{k)j, then n'{k) is statistically equivalent to n{k) with 
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E[n'2(A;)] = ct^ and E[n'j'^{k)] = Therefore 

y(k)=e^'-''/^''{l + n'j,{k)+n'j{k)j) 



= ^(1 + n'j^{k)y + n;2(A:)e(^'2"«/^''+^^(^)) (32) 
Where 9e{k) is phase noise corresponding to n{k). At high SNR, the following approximation holds: 

(33) 



9e{k) ^ tan{e,{k)) = -^^^^ « n'j{k) 



E 



l + n'j,{k) 



(34) 



Define y = [y(l), y(2), • • • y{N)] and ^ = [^{l), (^(2), • • • ip[N)], ^{k) = 27rq/Xk + 9e{k). It is clear 
that the original problem of estimating q using tp, corrupted by noise of variance cj^/2, is equivalent to 
estimating it from y with double variance. The Cramer-Rao bound of the equivalent problem is easily 
obtained as follows ll20l . 

Let Ofc = Iie{y{k)}, bk = Im{y{k)}. The probability distribution function is 

2N 



1 



exp 



^ iuk- cos{2'Kq/\k) +[hk- sin(27rg/Afc) 



CT; 



2N 



exp 



N 



— ^\[o-k- cos{2iTqfk/c) ) + [bk- sm{2TTqfk/c 



(J: 



k=l 



(35) 



The entry of Fisher information matrix and the Cramer-Rao bound (CRB) are given by 

'dlogf{y,q) 91og/(y,g)" 



E 



dq 



dq 



k 



k=l 



^2^2 / N 



CRB(g) = ([FU-i = ^(5:/, 



.2^2 / N 



Note that the CRB coincides with the HMSE, see (120] 



k=l 



47r2 



(36) 



(37) 



III. Frequency Design 

There exist three different frequency design goals, including maximizing UMR, minimizing the outlier 
probability and enhancing the MSE performance. We provide corresponding design criterion indepen- 
dently and then present a simple algorithm with all factors considered together. 
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A. Frequency Spacing Design for UMR and Outlier 

Proposition 2: Suppose the frequency resolution of the system is A/, that is to say, each frequency 
spacing must be integer number of A/. If all the frequency spacing can be expressed as A/j = PiAf, 
Pi is prime number, then the UMR is AL = c/A/. 

Proof: It is the corollary of Theorem [T] ■ 
Remark: Proposition |2] demonstrates that we can extend the measurement range by simply reducing the 
GCD and rearranging the elements of Af has no influence on the UMR. 

Define Ap = 27rAgAf^r^/c and Ap' = 2'n:AqAi'^/c . For a specific Af, assume the 



2lT 



2tt 



optimization problem in dTSl ) attains its minimal at Aq. From ( fTSl ). we can infer that Ap will then 
contain a large amount of value near zero. It is clear that the value of objective function will decrease 
as the number of extremely small value in Ap increases. Then a sharp peak, comparable with the main 
peak at the true qq, will create and the outlier probability increases. We then try to avoid the occurrence 
of extremely small value in Ap simultaneously. To this end, we provide the following proposition. 

Proposition 3: Assume the GCD of the frequency spacing Af is Af, i.e., A/j = kiAf, ki is positive 
integer. For any two frequency spacing A/j, Afj, if ki and kj are co-prime with GCD(/cj, kj) = 1, then 
for any given Aq, there exist at most \N + 1/2] elements in Ap which attain to zero, where [.] denotes 
the integer ceiling operation. 

Proof: If more than \N + 1/2] elements approaching zero, then there exist at least two pairs of 
consecutive zeros in Ap. Since 

Ap(A: - 1) = [27r(/fc - fi)Aq/c]2n 
Ap{k) = [27r(A+i - fi)Aq/c]2^ 
Ap'{k) = Ap{k) - Ap{k - 1) 
this suggests that the vector Ap' has at least two zeros entries, i.e.,Ap'(i) = Ap'(j) = 0. That is 

2TTAqAfi/c = iTTUi ^Aq = nic/{kiAf) 
2TrAqA fj/c = iTrnj ^ Aq = njc/{kjAf) 

Since GCD(A;j, kj) = 1, then Ui must have the factor ki. In other word, rii = rriiki, mj > 1 is positive 
integer. Similarly, rij = rrijkj. On the other hand, Aq < AL = c/A/, we get 
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Contradiction. ■ 

B. Frequency Spacing Design for MSB 

For fi S> B, which is typical in most measurement system, the HMSE expressed in (l20l ) is not 
sensitive to different frequency pattern. In other words, different frequency design methods have almost no 
(negligible) influence on the HMSE performance. Therefore, we will be concerned with the optimization 
of MMSE rather than HMSE, in the form of 

max(Af)^r^ (In-i - uu^/N^ TAf (38) 

Proposition 4: Suppose B2 > Bi, for any given measurement frequency in frequency band Bi with 
spacing Af = [A/i, A/2, • • • A/tv-i]'^, If we pick a new set of frequencies from band B2 in such a 
manner that Af = [A/(, A/2, • • • A/]y_^]^, A/^ > A/^, then the new measurement frequencies will 
result in high ranging accuracy. 
Proof: Let 

Af' = [A/{,A/^,...A/;,_i]^ 

= [A/i, A/2, • • • A/;v-i]^ + [di, d2, • • • cLn-i] 
= Af + d 

Where d = [di\i = 1, 2, • • • iV - 1, > 0]. Define Q = r'^ (iVIjv-i - uu^) T, then 

{N-i)j i>j 

Q{i,j) = 1{N i = j,Q{i,j)>0 (39) 

^{N i<j 

It follows that 

(Af')^r^ (iVIjv-i - uu^) TAf 

=(Af)^r^ (a^Itv-i - uu^) TAf + d^r^ (iVl7v-i - uu"^) Td 



>(Af)^r^ (iVI^v-i - uu^) TAf 
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Remark: For equal-spaced measurement frequency such as RIPS, let the frequency spacing be Af 
[A/, Af, ■ ■ ■ Af]^, it is easily verified that 



2 2 



\2-l C (Tfl 



47r2 Af^r^ (liv-i - uu^/N) TAf 
_ c^Uaj 
~ ATT^ApN{N^ - 1) 
_ cH2a^g{N-l) 
~ Air'^B'^NiN + 1) 

This means that the ranging accuracy of measurement method using equal-spaced frequency increases in 
proportion to as well as B^. 

Lemma 1: Let a = {a^, l<A;<A^|0<ai<---< a^} be a set of positive numbers sorted in 
ascending order and g = 7r(a) = [^fc, 1 < A; < A^] be a permutation of a. Denote the set of all the 
permutations of a as ^ and define the partial sums sequence b of g as 

k ^ N 

hk = Y.9^^<k<N,h = -Y,bk 

i=l k=l 

1 ^ 

k=l 

k=i \i=i k=i 1=1 / 
= /(g) 

where b and ^(b) are the mean and variance of b. The optimal permutation is defined as g*, satisfying 

1 ^ _ 2 

g* = argmax/(g) = argmax— V (bk - b) 
Then the optimal sequences are 

[ai, 03, as, . . . unjCln-i, ... 04, 02] N is odd 
[ai, as, as, . . . un-i, qn, ■ ■ ■ 0,4, 02] N is even 



g 



and 

[ai, a2, a4, . . . a^-i, cin, ■ ■ ■ 0,5, 03] is odd 



g 



[ai, a2, a4, . . . unjCln-i, • • • 05, 03] N is even 



Proof: We proved it in |[22l . 
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Lemma 2: Suppose the symbols a, g, g*, /(g) are defined as lemma. [T] Let g(^) = [fiffc+i, 1 < A; < 
— 1] and = {a^+i, \ <k < N — \] = {02, 03 • • • oat} denote the one-bit left-shift sequence of g 
and the subset of a. Define all the permutations of a^ as V = {v |v = 7r(as)}, then we have 

1) For any g and its left-shift sequence g(^), 

/(g) = (g«)'^r^ (i;v-i - uu^/iv) rgW 

2) The optimal permutation v* that maximizing (v)^r^ (^Iat-i — uu-^/A^^ Tv is just the one -bit 



left-shift sequence of g*, 



V* = arg max(v)^r^ (Iat.i - uu^/iV) Tv 

= M+i,i<fc<iv-i} = (g*f^ 

[03, as, . . . Oat, aAr_i, . . . 04, 02] is odd 
[as, as, . . . flAT-i, Cat, . . . 04, 02] is even 



and 



V = < 



[02, 04, • • • CLN-1: O-N: ■ ■ ■ 0^5) 0,3] is odd 

[02, 04, . . . Oat, OAf-i) • • • 05) 03] N is even 
Proof: 

1) For any g and its partial sums sequence b, let d = [di,d2,-" jC^Af]^ with dk 

J_ Y^A^ rl — X 

N ^k=l "fe — Af 



bk - gi, and 



d = ^ T,k=i dk = jf T,k=2 dk, then 



N 



1 



k=i \ ^ 

= E (^fc - (^91 + 1)02 + ■■■ gN)/N^ 

- N _ 2 _ " 



Define d(i) = [d2,d3,--- ,dN]'^ and g^^) = [c,2, c/3, • • • ,57v]'^, we can get 

u^rg(i) 



d(i) = rg«, d 



N 
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(d(^) - du) = rg 



(1) u^rga)u 



uu 



t-N-l 



N 



■ N _ , 



.k=2 



dW-duV fd«-du 



(g{i))^r^ (ijv-i - mm^/n) {1n-i - uu^ 



rg 



(1) 



/(g) = (g(i)) - uu^/iv) - uu^/n) + uu^/N^) rg(i) 

= (g«)'^r^(iA^-i-uu^/iv)rg« 

2) According to lemma [H gl = oi, then (g*)*^^'' G V, we have 

/(g*) = (g*('^)'^r^ (itv-i - uu^/iv) rg*(i) 
< (v*)^r^ (ijv-i - uu^/Af) rv* 

Let V* = {ai, V*} = |ai,i;i, • • • v* € V, therefore 

(v*)^r^ (l;v-i - uu^/iv) rv* = /(v*) < /(g*) 
From (l40l) and (|4T]) . it is readily seen that 

(v*)^r^ - uu^/iv) rv* = /(g*) 
V* = (g*)(^) 



(40) 



(41) 



Based on lemma ID we deduce the following Theorem 

Theorem 2: Suppose that Af = [A/i, A/2, • • • A/jv_i] is any given frequency spacing sorted in 
ascending order, then the optimal permutation of Af, in the sense of maximizing equation(l38l). takes 
the form 

Af* = [A/i, A/3, • • • A/^_i, A/;v-2, • • • A/4, A/2] 

Proof: It is a direct conclusion of lemma |2] ■ 
We named the permutation of the form Af * provided in Theorem |2] as min-error permutation and the 
dual form Af = [A/tv-i, • • • A/3, A/i, A/2, A/4, • • • A/7V-2] as max-error permutation. It is proved 
in |[T5l that the permutation in the form of Af has extremely small variance once substituted into 
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equation(l38l). This implies that this permutation will lead to large ranging error when used in measurement 
frequency design, (use {} to denote set) 

Theorem 3: For a measurement system with bandwidth B, frequency number N, frequency resolution 
A/mm> the total frequency number [M = B/Af^in\ i[\ denotes the floor function), if the range fall in 
[—c/2B,c/2B], then the optimal measurement frequencies must be picked up from as near as possible 
to both ends of the frequency band, with frequency spacing in the form of 

Af* = [1, 1, • • • M + 2 - iV, • • • 1, ifAU,^ 

Proof: See Appendix iBl ■ 

C. Algorithm 

Based on Proposition |2H1] and Theorem |2l we proposed a quite simple yet effective(modify) algorithm 
in Table H 

Remark: The large enough primes set in Table U allow the algorithm to select the appropriate primes as 
needed. In the latter part of the paper, the notation {B, N, A/min, i, K) stands for the design parameters 
of the min-error method defined above. Note that the UMR of the method is AL ^ c/{KlS.f^\^). The 
common factor K and the index i have great influence on the UMR as well as the estimation accuracy, 
as is shown in the following simulation. 

IV. Simulation Results 

In this section, we present simulation results to compare different frequency design methods under 
various scenarios. The measurement frequencies of the Towers method satisfy Bl 

/s ~ fl _ fi — fl _ _ fl _ h 

/2 - fl /s - fl fN - fl B 

The measurement frequencies of RIPS are IJl 

fi = fi + {i-l)B/{N -l),i = l,2---N 
The frequencies of the constrained optimal method are given by 

Af = [1, 1, . . . , M + 2 - iV, . . . , 1, ifAUin 

i-l 

fi = fi + Yl Af (fe), M = i?/A/,,i„ 

k=l 
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* Define the set of primes less than M as Vm- The positive integer 
M is selected so that IVm] > iV. 

* Assume A/min is the frequency resolution of measurement system 
and K is the common factor, i.e., KAf^in is the minimal frequency 
spacing of the frequency. To make full use of the bandwidth B, we 
perform the following steps (assimie the UMR requirement is known 
as priori information ) : 

1) i = l. 

2) Pick the iV — 1 consecutive primes subset 5' = 
{VM{i) ■ ■ -VMiN + i — 2)} and find the correspond- 
ing K, obey KAUi^Y^jJi ^'U) < B and (K + 

3) if c/(/f A/min) > UMR, go to step 4); else, i = i -|- 1 and go 

to step 2). 

4) Then, the frequency spacing set AJ^ — {Afj — 
KAU^S^j) |j = 1 • • • TV - 1} satisfying J2jJ,' ^fi < B 
is constructed. 

■k Sort the frequency spacing in ascending order and obtain Af = 
[A/i,A/2,...A/;v-i]. 

• If A'^ is odd, Af can be rearranged as 

Af* = [A/l, A/3, • • • AfM-2, AfM-i,- ■ ■ A/4, A/2] 

• If iV is even, then 

Af* = [A/i, A/3, • • • A/;v-l, A/;v-2, • • • A/4, A/2] 

* The measurement frequencies is finally obtained fi = /i + 
E;-\Af*(fc),i = l,2,...iV. 



The frequencies of the prime-based min-error method are selected as (We assume N to be even) 

Af = [A/i, A/3, • • • A/;v-i, AfN-2, ■ ■ ■ A/4, A/2] 
fi = h + Y,^f{k) 

k=l 
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To illustrate the effect of permutation, we also provide the worst permutation used in the prime-based 
max-error method(with the same set of frequency spacing but a different permutation), in the form of 

Af = [A/jv-i, A/jv-3, • • • A/i, A/2, • • • AfN-4, A/iv-2] 

j-i 

/i = /i + ^Af(A;) 

k=l 

Unless otherwise mentioned, the following parameter values are assumed: The measurement frequencies 
range from 400 MHz to 500 MHz, the initial frequency fi = 400 MHz is used. The interferometric range 
is qo = Om, the phase error in each frequency is modeled as independent and identically distributed 
(i.i.d.) zero-mean complex white Gaussian noise with variance E{6'^{k)} = Uq. The signal-to-noise ratio 
is defined as SNR = l/(2cr|) according to section UPEl We adopt the LS search algorithm to find the 
optimal solution, with a step size of 0.001 m ranging from — c/2A/ to c/2A/, where Af = B/{N — 1) 
and c/A/ are the frequency separation and UMR of RIPS method, respectively. 

The MSE of different methods are averaged over 2000 Monte Carlo runs for each SNR. The MMSE 
of RIPS method. Prime-based min-error method and Prime-based max-error method are denoted by 
MMSE-RIPS, MMSE-min-error and MMSE-max-error. The CRB of all the three methods are nearly not 
distinguishable, only the CRB bound of RIPS method is plotted using (|37] ) and denoted as CRB-RIPS. 

For a fair comparison among different frequency design methods, the same number of frequencies 
N and bandwidth B have been used. The notation {B,N,/S.f^\^,i,K) is defined in section [TlI-CI with 
A/min = 65 Hz is assumed, which is the frequency resolution of MICA2 platform used in the field 
experiment of section |V][r|. With no information about UMR requirement, i = 1 is used in the simulation 
except FiglT] Note that even in this case, the UMR is AL c/(A'A/min) > [c/ B)Y,^^i S'^{j) > 
{c/B){N - 1) = ALrips with ALrips denotes the UMR of RIPS, see Table H 

Figj3al Figj3b] illustrate that the constrained optimal method has the best ranging accuracy if the prior 
information that go is in the range [—c/2B, c/2B] is provided. It will fail to work once go is outside this 
region. This feature is predicted in Theorem [3] and reduces its measurement range greatly. The proposed 
min-error method outperforms all the others and is slightly inferior to the constrained optimal method 
only in the above limited measurement range. The MSE performance of towers method is not satisfactory 
due to the fact it uses a local instead of global method to optimize its frequency. 

When the condition /i S> -B is met, the multi-frequency interferometry ranging technique exhibits the 
unique double-threshold feature, as shown in FigJH That is to say, besides the classical threshold, another 
threshold occurs. The MSE curve firstly follows the MMSE derived in (|38] ) tightly once SNR exceeds 
the classical threshold and drops to the CRB (or HMSE) when SNR reaching the second threshold. Only 
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SNR(clB) 



Fig. 4. The double MSB curve for = 21 and B = 20 MHz. 



the CRB of RIPS is plotted for minor difference between the three methods and that is just the reason 
of optimizing MMSE in the paper. The prime-based min-error method has the minimal MMSE as well 
as the best actual accuracy as expected. Compared with the prime-based max-error method, the merit of 
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Fig. 5. The impact of N and B. 



frequency spacing rearrangement is clearly visible, with about 4 dB gain obtained at the MSE of 10~ . 
This point will be further illustrated in later simulation. 

From Fig|5l it is observed that increasing the number of frequency or the measurement bandwidth 
will both improve the ranging accuracy of all the methods. The performance improvement achieved by 
increasing bandwidth B is more significant, relative to the increment of N. Note that the proposed prime- 
based min-error method also works well for relatively wide bandwidth, also seen in Figj3l although the 
design principle is derived under the assumption of fi » B. 

With fixed B and N in Figj6j it is interesting to find that the CRB is more easily attained at the cost 
of larger CRB for a low initial frequency fi. The opposite has be seen for a high initial frequency. This 
phenomenon coincides with the theory. 

The influence of the parameter i,K on the MSE performance for a particular SNR is shown in FigjTJ 
with = 31 and B = 40 MHz. The x-coordinate is the combination of i,K, which is the abbreviation 
of (i?, N, A/min, i, K) for simplicity(i is always set to 1 except in this case). The superiority of the min- 
error permutation method is clearly seen for large K and will diminish as K decreases. The reason is 
that the difference between the frequency separation of the proposed method decreases and the frequency 
pattern tends to approach the RIPS method with decreasing K. 

When SNR is below the conventional threshold and the outlier occurs, the error is uniformly distributed 
across the entire measurement range and the MSE will not reflect the estimation performance any 
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10 



10 



N=21,B=20MHz 




-^1^ — Prime-based min-error method 
-a — RIPS method 
CRB-RIPS 



10 15 20 

SNR(dB) 



25 



30 



Fig. 6. The impact of tlie initial measurement frequency. 



10 



SNR=1dB 




SNR=5dB 



— N=31, B=40MHz, prime-based min-error method 
-B — N=31, B=40MHz, RIPS method 



10"' 

(1,386) (2,358) (4,311) (8,243) (16,166) (32,98) (64,51) (128,25) 
impact of (i,K) 



Fig. 7. The impact of the common factor K and the prime index i in algorithm IIII-CI 



more(one large outlier may ruin the MSE curve and lead to larger MSE than that caused by many small 
outliers, so the prime-based min-error method may be inferior to RIPS method in MSE performance for 
low SNR, see FigJJ]). Hence, we compare the probability of incorrect unwrapping of different methods 
rather than MSE in FigH] From (dJ and dD, the incorrect unwrapping probability Pj is defined as (for 
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SNR(dB) 



Fig. 8. The probability of incorrect phase unwrapping. 



RIPS method jhe min-error method 




q (m) q (m) 



(a) (b) 
Fig. 9. The unambiguous measurement range versus frequency pattern, qo = lOOm.(a) RIPS method.(b) the min-error method. 



correct unwrapping, the error must be less than one wavelength) 

Pf = P(n 7^ n),n = [ni,n2,- • • ,nAr],n = [ni,n2,- • • ,nN] 
hi = round(g/Ai — (p{i)/2Tr),ni = round(go/Ai) 
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where round(x) returns the nearest integer to x and q is the estimation of go, the results is averaged 
over 500,000 Monte Carlo runs. FigH] shows that the proposed min-error method has a much better 
performance, in term of Py, than both the RIPS method and the max-error method. The min-error method 
provides approximately 2 dB and 4 dB gain over the RIPS method and the max-error method respectively, 
at an incorrect unwrapping probability of 10~^ for N = 21,B = 20 MHz. Similar observation also holds 
for AT = 41,5 = 40 MHz. 

Figj9] further compares the UMR of the min-error method with RIPS method for the above parameter 
N = Al,B = 40 MHz. The design parameter of the proposed method is (40,41,65,1,199). So the 
UMR of RIPS method and the proposed method are AL = c/{B/{N - 1)) = 300m and AL = 
c/(i^A/min) — £Ylif=i I Si^i ~ 23193m, which both are in good agreement with simulation 
results. The proposed method achieves far more large UMR than RIPS. From Fig|9bl it is also noted that 
larger sidelobes are not seen for being uniformly averaged over the whole parameter space and sharp 
peak is observed at the true location. These properties result in superior estimation accuracy and low 
outlier probability verified by the MSB curve. Moreover, the UMR can be easily enlarged by adjusting 
the parameter K. 

V. Experimental Evaluation 

The low-cost mica2 nodes is exploited for field experiment. All the ranging procedure is similar to 
the one in |1| and |2| except the measurement frequencies. Five nodes are used with two transmitters 
(A and B) and three receivers (C~E). For each measurement round, two transmitters and two receivers 
are needed. So the nodes form three deployment scenarios, ie.ABCD, ABCE and ABDE. The nodes are 
deployed in football-field of our campus, see Fig[TOl To alleviate the multipath effect, all the nodes are 
placed one-meter above the ground, seen ||2l| for details. The real coordinates of nodes are determined 
via differential GPS, which has position error of about ±2 cm. The true ranges are (Iabcd = 19.19m, 
dABCE = 6.88m and cIabde = —12.31m, where cLabcd is the linear combination of the distances in 
the form of (Iabcd = dAD — Abd + ^bc — dAc< dxY denotes the distance between node X and Y. 

The measurement frequencies range from 410MHz to 450.378MHz, the bandwidth and frequency 
number are B = 40.378 MHz, = 31. The following methods are compared under identical B and N: 
RIPS, min-error method, max-error method and the random method. For random method, the measurement 
frequencies are picked randomly from the usable frequency band. The experiment parameter of prime- 
based min-error method is {B, N, A/min, i, K) = (40.378, 31, 65, 12, 200). Since i = 12, then the prime 
sequence described in section UlTC] is 5 = [37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 
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.02) 










(23.85,52.53) ^^^^ (42.15,52.50) 




(53.15,16.50) 

•d 



(65.99,0) 



Fig. 10. Coordinates of nodes (m). 



113 127 131 137 139 149 151 157 163 167 173 179]. The nodes have fine frequency resolution of A/mm = 
65 Hz. It is easy to verify that J2i ^ (O-^^/min = B holds. Based on Propositon [TJ the practical UMR 
is AL = c/ (KA/min) - eEili ^ Eili ~ 23.077km. For each method and deployment, 50 
independent experiments are performed. 

Since the UMR is AL ^ 23.077km, the search ranging is set to [-1000m,24000m] to guarantee one 
ambiguity solution being included in Fig[TT] The error between the real and estimated range is then 
plotted in FigfTTl 

It is obvious that the errors lie in two regions, one is near zero and the other is located at 23077 m or 
23078 m. The latter is very close to the practical UMR(P-UMR). This implies that the P-UMR holds true 
and is robust to noise and frequency inaccuracy, which are inevitable especially for low-cost hardware 
such as mica2 node. 

The search operation of the following experiment is performed in the UMR range of RIPS, ie. 
[— c/2A/, c/2A/] [—100 m, 100 m] to avoid ambiguity. The distribution of ranging error is shown in 
Fig[T2] for different methods. Fig{T3] exhibits the cumulative distribution function (CDF) of the absolute 
value of error for the deployment of ABDE. These two figures confirm the results in FiglS] It is pointed 
that the ranging errors are bias and not Gaussian distribution any more due to the existence of multipath 
or wrapped Gaussian noise. Even in this case, the proposed min-eiTor method is still superior to the other 
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Fig. 11. Validation of the unambiguous measurement range via field experiment, (a) the min-error method, (b) the max-error 
method. 



schemes. 

The MSE performance of different deployments are also shown in Fig{l4l The discrepancies are 
observed for all the methods except the min-error method, which has the best accuracy all the time. This 
discrepancies can be attributed to the multipath effect. It is well known that frequency-selective shading 
is introduced in multipath channel and it varies from one location to another. This is so-called frequency- 
selective and space-selective property of wireless channel. Therefore, the measurement frequencies of a 
certain method may undergo deep shading in one deployment and leading to small SNR in the receiver. 
It is quite possible that the opposite happens for another deployment. Consequently, the discrepancies 
appear. 

VI. Conclusion and Discussion 

In this paper, we focused on the frequency optimization of MFI system to extend the UMR, decrease 
the probability of outlier and improve the estimation accuracy. We prove that the UMR of MFI is in 
inverse proportion to the greatest common divisor (GCD) of frequency interval, thus reducing the GCD 
will extend the measurement range. We also point that the theoretical UMR of MFI is usually optimistic 
and provide an expression for the practical UMR. We explore the relationship between outlier probability 



February 8, 2012 



DRAFT 



SUBMITTED TO IEEE TRANS. X X, VOL. X, NO. X, X 2012 



30 



deployment: ABDE 



deployment: ABDE 




0.2 0.4 0.6 
error (m) 

(a) 

deployment: ABDE 



(b) 

deployment: ABDE 




-0.5 
error (m) 

(c) 



0.2 0.4 0.6 
error (m) 



(d) 



Fig. 12. Error distribution of range estimation with 50 independent experiments for each method, (a) min-error method, (b) 
max-error method, (c) RIPS method, (d) random method. 



and the ambiguity function and suggest to use prime-based frequency interval for outlier suppression. 
The unique "double threshold" phenomenon of MFI using LS estimator is discovered and the expressions 
for the MMSE and HMSE are derived. Focusing on the optimization of MMSE performance, we provide 
an optimal permutation for any set of frequency interval and prove its optimality. Based on the finding 
mentioned above, we present a quite simple and effective frequency design method. Simulation results 
verified that the proposed method outperforms the existing method in UMR and MSE simultaneously. 
Field experiments further demonstrate the new method's robustness to practical interference such as 
frequency inaccuracy and multipath error. Although the theory is developed under the assumption of 
narrow-band ranging with fi S> B, it is found that the proposed method also performs well in relatively 
wide bandwidth. As a final remark, since the outlier suppression by non-searching-based frequency 
interval design is an open problem in the literature, we just give a heuristic design method and further 
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Fig. 13. CDF of the absolute errors for different methods. 
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Fig. 14. MSE for different methods with three deployments, averaged over 50 independent experiments for each method and 
deployment. 

research is needed in this area. 
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Appendix A 
Proof of the Proposition [H 

Suppose the true range is qq. Since /j = {ki + e)A/min, it is obvious that the cost function could 
not achieve zero at qq + AL, AL = c/A/mm- We will instead consider the location of local minimal 
Q = Qo + AL + X, near the original ambiguous location. The cost function is rewritten as 

N 



i=l 
N 

= E 



i=l 

N 



2^So _ ^^(go + AL + x) 
A,; A,: 



2n 



E 

i=l 

N 

E 

i=l 

T{x) 



27rx , 2^((fci + e)A/mi„) AL 



A,: 



+ 



27r 



27r(x + eAj 



A,; 



2ix 



(42) 



we notice that T{x) will decrease when x has opposite sign relative to eAj. Assume Ai > A2 • • • > Aat 
and Ai < 2Aj (it is reasonable considering fi » B), If < |e|Ai, then \x + eAj| < Aj/2. The problem 
is to find the local minimal x within the constraint Ixl < lelAi, we have 



^ /2^(x + eAi'^ ^ 



Therefore, 



Then 
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We obtain 



N 



T{x) < ^ (27re(Ai - A;v)/Atv)' = ^Ntt^'' {h/B) 



(46) 



(47) 



i=l 



February 8, 2012 



DRAFT 



SUBMITTED TO IEEE TRANS. X X, VOL. X, NO. X, X 2012 



33 



where B = — fi. we, finally obtain 



lim S{q\q = qo + AL + x) = 

/i/B->oo 

lim S{q\q = qo + AL + x) = 
e— >0 



Appendix B 
Proof of the Theorem [3] 

We firstly sort the — 1 adjacent frequency spacing in ascending order 

Af = [A/i, A/2, • • • AfN-i] = [h,k2,--- kN-i]AU 

= M 



(48) 



N-l 

Y.h = B/AU 

i=l 



l<ki<k2<--- kN-i,kN~i <M + 2-N 

According to Theorem |2l the optimal rearrangement of the spacing is 

Af = [A/i, A/3, • • • A/jv-i, A/;v-2, • • • A/4, A/2] 

= [h,k3, ■ ■ ■ kN-i,kN-2, ■■■ki, /i;2]A/mm 

If kf^_i < M + 2 — N , then there exists at least one k,n> 2,m ^ N — 1. Without loss of generality, N 
and m are assumed to be even and odd). Then 



Af 
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Note that Q{i,i) + Q(j,i) > 2Q(i, j), then 

'm + 1 m + 1 



tr (AuAu^Q) = Q ( ) + Q 

> 

Let U = AuAf^, since Q is symmetric, we have 
tr (Af Au^q) = tr (AuAf^Q) 
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By Theorem H Af(j) < Af(iV - 1 - j), therefore 
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N/2-1 



m + 1 



m + 1 .' 



+ ^ Af(i) jN 
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i=(m+2)/2 



>0 



> 



It follows that 



Af'^QAf < Af'^QAf' 



(49) 



This implies that the sequence Af is superior to Af. Therefore, the maximum value kjq^i of the 
optimal frequency spacing Af*, normalized by A/min, must satisfy /cat-i = M + 2 — N . Otherwise, for 
any Af with kj^_i < M + 2 — N, we can always find another frequency spacing better than it on the 
basis of ( |49l) . 

The key idea of the proof is that all the Af may be classified into different categories (categorized 
into groups). Those Af composed of the same set of frequency spacing belong to the same class, which 
differ only in the permutation. For each class, the optimal one is easily obtained by Theorem |2] The 
global optimal across all the local optimal takes the form: 



It is pointed that when the frequencies are densely distributed on both ends of the band, those 
frequencies can be approximated as two frequencies on both ends and the corresponding UMR approaches 
c/B. Thus, the condition of g € [—c/2B,c/2B] is imposed on the Theorem to exclude ambiguity. 
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